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Abstract 

Explicit Runge-Kutta schemes with large stable step sizes are developed for inte- 
gration of high order spectral difference spatial discretization on quadrilateral grids. 
The new schemes permit an effective time step that is substantially larger than the 
maximum admissible time step of standard explicit Runge-Kutta schemes available 
in literature. Furthermore, they have a small principal error norm and admit a low- 
storage implementation. The advantages of the new schemes are demonstrated through 
application to the Euler equations and the linearized Euler equations. 

1 Introduction 

Throughout the past two decades, the development of high-order accurate spatial discretiza- 
tion has been one of the major fields of research in numerical analysis, computational fluid 
dynamics (CFD), computational aeroacoustics, computational electromagnetism and in gen- 
eral computational physics characterized by linear and nonlinear wave propagation phenom- 
ena. High-order discretizations have the potential to improve the computational efficiency 
required to achieve a desired error level by allowing use of coarser grids. Indeed, in modern 
wave propagation problems characterized by complicated geometries, complex physics and 
a wide disparity of length scales (e.g. large eddy simulation, turbulent combustion, flow 
around flapping wings, rotor-blade interaction), the need for high-accuracy solutions leads 
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to a prohibitive computational cost when low-order (i.e. first- and second-order accurate) 
discretizations are used. High-order schemes have much better wave propagation proper- 
ties and a truncation error that decreases more rapidly than that of low-order schemes if the 
solution is sufficiently smooth. Therefore, for problems that require very low numerical dissi- 
pation and small error levels, it may be advantageous to use high-order spatial discretization 
schemes; see for instance [H El El EH]- 

Among high-order methods, the spectral difference (SD) scheme is receiving increasing 
attention P [131 EH ES ED EHl ES]. The SD scheme offers several interesting properties. It 
is capable to obtain solutions with arbitrarily high order of accuracy. It can be applied to 
unstructured quadrilateral and hexahedral meshes (tensorial cells). The conservation laws 
to be solved are in differential form, making implementation of new equations easier, and 
avoiding the use of costly high-order accurate quadrature formulas. 

Although the formulation of high-order spatial discretization is now fairly mature, the 
development of techniques for efficiently solving systems of ordinary differential equations 
(ODEs) arising from high-order accurate spatial discretizations has received less attention. 
The cost of solving an initial value problem up to a fixed time is inversely proportional to 
the time step used, so it is desirable to use the largest step size possible if the temporal 
discretization error is acceptable. For higher order schemes, the spectrum of the Jacobian 
of the semi-discretization often has increasingly large eigenvalues. As a result, the step size 
is often limited by stability requirements, which become stricter with higher order methods. 
Implicit methods allow the use of much larger step sizes, but lead to very large memory 
requirements that may not be feasible. By using explicit methods with large number of 
stages, such memory requirements can be avoided and larger stable step sizes can be achieved. 

This work focuses on the development of new optimized explicit Runge-Kutta (ERK) 
schemes to compute wave propagation efficiently and accurately with high-order SD meth- 
ods on unstructured uniform or quasi-uniform non-simplex cell grids. The schemes are 
optimized with respect to the spectrum of the SD discretization, using the two-dimensional 
(2D) advection equation as a model problem. Linear stability optimization determines the 
coefficients of the stability polynomial but does not uniquely determine the full RK method. 
A second optimization step is used to determine the Butcher coefficients of the scheme, op- 
timized for a small leading truncation error constant and low-storage form. The low-storage 
form is crucial for memory reasons, since many stages are used. 

Many authors have studied the design of optimal ERK schemes with many stages for 
integration of high order discretizations of partial differential equations (PDEs). Past efforts 
focused on schemes with a relatively smaller number of stages [21 ESI HH [291 E] • By using 
the algorithm developed in |25, we are able to develop schemes with much larger number 
of stages and with higher order of accuracy. Our work is also the first to develop schemes 
specifically for the spectral difference semi-discretization. Whereas past studies have focused 
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on application to linear problems only, and typically employed structured grids, we validate 
the effectiveness of our methods also on a nonlinear, fully unstructured example. Our new 
optimal ERK methods increase the computational efficiency of the SD method for wave 
propagation problems up to 65% and 57%, respectively for 4*'^-order and 5*'^-order accurate 
simulations. 

The remainder of the paper is organized as follows. In Section|2| we review the SD method 
for non-simplex cells. Section |3] is devoted to the description of the two-step optimization 
algorithm used to design new ERK schemes, in which we first select an optimal stability 
polynomial and then design a corresponding ERK method. We also present the main features 
of the optimized methods and discuss their theoretical efficiency. Section|4]presents numerical 
results for three benchmark test problems, which demonstrate that the new schemes lead 
to large performance gains over standard ERK schemes available in the literature, for both 
linear and nonlinear problems. Conclusions and future directions are given in Section |5| 



2 Spectral difference discretization 



In this section, we review the spectral difference approach to semi-discretization of hyperbolic 
conservation laws. 

Consider the general hyperbolic system of conservation laws over a rf-dimensional domain 
C M"' with boundary dQ and completed with consistent initial and boundary conditions: 



^ + V.f(q) = 




s(qj 



on Q 
on dQ, 



Here, x, q, f, s, t*^ and are respectively the position vector, the vector of the conserved 
variables, the flux vector, source terms, and the lower and upper bound of the time interval. 
The spatial domain Q is discretized into tensor product cells (quadrilateral and hexahedral 
cells) with domain and boundary Qi and dQi. 

For each cell i, take a mapped coordinate system = r], Q . The transformation from 
the standard to the physical element in the global Cartesian coordinates for the cell i is given 
by 

Vi (C, V, C) 



Xi 



(2) 



with Jacobian matrix jj and Jacobian determinant Jj. The fluxes projected in the mapped 
coordinate system (f/) are then related to the flux components in the global coordinate 
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system by 



















. hi . 



7 f 



(3) 



Therefore, the hyperbohc system ([T]) can be written in the mapped coordinate system as 



9t 



dt 



dr] d( 



(4) 



where q/ = Jjq and V ^ are the conserved variables and the differential operator in the 
mapped coordinate system, respectively. 

For a (p + l)-th-order accurate d-dimensional scheme, solution collocation points with 
index j are introduced at positions ^| in each cell i, with A^^ given by N' = {p + If. Given 
the values at these points, a polynomial approximation of degree p of the solution in cell i can 
be constructed. This polynomial is called the solution polynomial and is usually composed 

of a set of Lagrangian basis polynomial Lj (^] of degree p: 



l,...,iv^ 



(5) 



Therefore, the interpolation coefficients are given as Qij = Qi y^jj where Qij are the 
conserved variables at the solution points, i.e. the unknowns of the SD method. 

The divergence of the mapped fluxes V ^ ■ f ^ at the solution points is computed by in- 
troducing a set of N-^ flux collocation points with index / and at positions , supporting 

a polynomial of degree p + 1. The evolution of the mapped flux vector f ^ in cell i is then 
approximated by a flux polynomial f/, which is obtained by reconstructing the solution vari- 
ables at the flux points and evaluating the fluxes f/^ at these points. The flux is represented 
by a separate Lagrange polynomial: 



Sir 



Lm = 1 N^. 



(6) 
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Hence, the coefficients of the flux interpolation are defined as 



where F^^^ is the flux vector at the cell interface. In order to maintain conservation at a cell 

level, the flux component normal to a face (i.e. F^^^^ ■ n^) must be continuous between two 
neighboring cells. However, the solution at a face is in general not continuous and requires 
the solution of a Riemann problem. For many nonlinear hyperbolic systems, such as the 
compressible Euler equations, the exact Riemann solution cannot be written in closed form 
and is prohibitively expensive to compute. Therefore, cheaper approximate Riemann solvers 

are typically used. The tangential component of F^^^ is usually taken from the interior cell 
(see for instance |12]). 

Taking the divergence of the flux polynomial V ^ ■ F^ in the solution points results in the 
following modified form of (|4]), describing the evolution of the conservative variables in the 
solution points: 



dt 



V-F,; 



v«-f/ 



where Fj is the flux polynomial vector in the physical space whereas R^j is the SD residual 
associated with Qij- This is a system of ODEs, in time, for the unknowns Qi.j. 



2.1 Solution and flux points distributions 

Huynh (TH] showed that for quadrilateral and hexahedral cells, tensor product flux point 
distributions based on a one- dimensional (ID) flux point distribution consisting of the end 
points and the Legendre-Gauss quadrature points lead to stable schemes for arbitrary order 
of accuracy. 

In 2008, Van den Abeele et al. showed an interesting property of the SD method, 
namely that it is independent of the positions of its solution points in most general circum- 
stances, for both simplex and tensor-product cells. In the above work it has been shown that 
the distribution of the solution points has very little influence on the properties of the SD 
schemes, and in fact, for linear problems, different distributions lead to identical results. This 
property greatly simplifies the design of SD schemes, since only the flux point distributions 
has to be taken care with. It also implies an important improvement in efficiency, since the 
solution points can be placed at flux point positions and thus a significant number of solution 
reconstructions can be avoided. Recently, this property has been proved by Jameson [2U]. 
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Figure [T] shows a typical distribution of flux and solution points for a third-order SD 
scheme in 2D. 



• ^ *- 



-1 



Figure 1: Typical 3'''^-order {p = 2) quadrilateral SD cells, with component- wise flux point distri- 
butions. Solution (o), ^- (▼) and ry-flux points (a). 



2.2 Advection equation 

The advection equation represents the simplest hyperbolic conservation law. It models the 
advection of a scalar conserved variable q with constant advection speed a. The conserved 
variables and the convective flux are then 

q = 9) 

. : (9) 

I = aq. 



Therefore, the conservation law reads 

dq 



+ V ■ {aq) = 0. (10) 



2.2.1 Spectrum of the two-dimensional advection equation 

Despite its simplicity, the one-dimensional advection equation is often used as a model prob- 
lem to design and analyze new spatial discretization schemes for convection-dominated prob- 
lems. It is also used to optimize the coefficients of time integration algorithms for CFD, and 
for wave propagation problems in general. 

In this work, however, we use the 2D advection equation as a model. In the 2D case, the 
discrete operator arising from the spatial discretization is also a function of the convective 
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velocity direction. This approach allows us to consider different wave propagation trajectories 
and optimize the RK coefficients by using a richer spectrum (or Fourier footprint) than that 
of the ID advection equation. The richer spectrum leads to a design of more robust schemes. 
In the remaining part of this section, the procedure used to compute the Fourier footprint 
is described. 

Equation (10) is discretized in space by the SD scheme. A uniform grid with periodic 
boundary conditions is considered. The grid is defined by a generating pattern, which is the 
smallest part from which the full grid can be reconstructed by periodically repeating the 
pattern in all directions. For the 2D case and uniform quadrilateral meshes, the generating 
pattern is completely defined by the vectors ri and r2 (see Figure |2]) whose non-dimensional 
form is obtained by scaling them with the length of fi, denoted by Ar: fi = Arf[ and 
r*2 = Arr2- If the dimensionless vector r[ is chosen as [1 0]"^, then the dimensionless mesh 
is completely defined by the two components of rg. 




Figure 2: Generating pattern. 



The advection speed a in Equation (10) is defined by its amplitude a and orientation 
angle ip: 

cosip 



a = a 



At cell faces the solution is discontinuous, so two values for the convected variables are 
available. The normal flux component is calculated using following approximate Riemann 
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solver: 



F{QL,QR)-ln = a-l 



Ql + Qr 



d ■ ij 



Qb. — Ql 



:i2i 



where !„ is the unit normal oriented from the left to the right side and indices L and R 
indicate respectively the left and right neighboring cell to a face. In the present analysis, 
the internal component in each cell is used for the tangential flux component. 



After the SD semi-discretization of (10) on a uniform quadrilateral mesh, the following 
system of ODEs is obtained: 



dt 



-iO,+l 



(13) 



where the five matrices T are determined by the coefficients of the spatial discretization. 
They depend on the order of accuracy p of the SD scheme, the generating pattern, and the 
advection velocity orientation angle ip. The column vector Qj j contains all solution point 
variables of the cell with indices i and j (Figure [2]). 
Inserting the following plane Fourier wave 



Q(t) e^H^i+^^'i) 



(14) 



into Equation (13) results in 

dq, 



dt Ar ^ 



(15) 



where / = is the imaginary unit number. Here, k and K are the wave vector and the 
dimensionless wave vector given by 



k = k 



cos 6* 
sin 6 



and 



K = kAr = K 



cos 6 
sin^ 



(16) 



(17) 
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Equation (15) can be written as 



where the matrix L is defined by the SD spatial operator. The set of eigenvalues of L is 
the spectrum or Fourier footprint of the spatial discretization. The spectrum depends on 
the order of accuracy p of the SD scheme, the generating pattern, the direction ip of the 
convective velocity, and the dimensionless wave number vector K. 
Here we take a uniform Cartesian grid defined by 

^ r-i^( a... ) (19) 



y ' ' ^ 

and, for a given order of accuracy, we compute the spectrum of the operator L by varying 
ip, K and 0. 



3 Optimized Runge— Kutta schemes 

The spectral difference semi-discretization of a PDE described in the previous section leads 
to an initial value problem 

' Q'(t) = F(Q) 



Q(o) = Qo ' 

where Q(t) : M — ?► ]R^°°^ and F : ]R^°°^ — )■ M^^^'' are the vector of the unknowns and 
the vector of the residuals, respectively. The number of degrees of freedom is denoted by 
jyDOF _ X A^*, where N is the number of cells used to discretize the domain VL. System 



(20) is typically integrated by using a high-order accurate ERK time discretization, which 
takes the form 

i-1 

= Q" + At^ai,F(r,) (21) 
i=i 

i-1 

gn+1 = gn + At 5^ hF{Yj), (22) 
i=i 

for a scalar ODE. The properties of the RK method are determine by its coefficient matrix 
A = [aij\ and column vector b = \b.j\ which are referred to as the Butcher coefficients [5] . In 
this section, we describe our approach to designing RK schemes that maximize the absolutely 
stable time step size, have reasonably small error constants, and can be implemented with 
low storage requirements. 
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3.1 Optimization of the stability polynomial 



Stability of RK integration is studied by applying the method (21) to the linear scalar test 
problem Q'{t) = X Q. Any RK method applied to this problem yields an iteration of the 
form 

g"+^ = ^{At A) g" (23) 

where the stability function ip^z) depends only on the coefficients of the RK method ( [T5| 
Section 4.3] [IS]): 

s 

^{z) = l + J2 b^A^'^^e z^. (24) 

3=0 

Here e is a column vector if size s made by ones. The stability function governs the local 
propagation of errors, since any perturbation to the solution will be multiplied by ^/'(At A) 
at each subsequent step. 



We say the iteration (23) is absolutely stable if 



At XeS where S = {z e C : \ip{z)\ < 1}. (25) 

The set S is referred to as the absolute stability region. 

When applied to a linear system of PDEs (such as the advection or linearized Euler 
equations discussed in the previous section), the SD semi-discretization leads to a linear. 



constant-coefficient initial value problem; i.e. (20) with F(Q) = ^ L Q where L is a 
fixed square matrix. Then application of a RK time integration again leads to the iteration 
(23), but with At A replaced by uh, where u = a ^ is the CFL number. Assume that L 



is diagonalizable and let Xi,i = 1, . . . , N^^^ denote its eigenvalues. Then the solution is 
absolutely stable for CFL number u if 

uXieS for 1 < i < iV°°^. (26) 

Thus the maximum absolutely stable step size is 

z/stab = max{z/ > : |^(z/Ai)| < 1 for ^ = 1, . . . ,iV°o^}. (27) 

Although this analysis is based on the linear problem, it is often used to obtain a practi- 



cal step size restriction for the nonlinear problem (20) by considering the spectrum of the 
Jacobian of F(Q). 
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In general, ip{z) for an s-stage, order p, ERK method is a polynomial of degree s that 
differs from the exponential function by terms of order z^^^: 

^{z) = j2 p^^' = E y + E (28) 



Comparing (28) with (24), we see that (5j = h'^ A^-^e. 

It is natural then to design optimal polynomials by choosing the coefficients Pj in (28) 
so as to maximize z/gtab- The optimization problem may be stated formally as follows 

Problem 1 (Stability polynomial optimization) 

maximize 
subject to 

|V'('^A)|<1 forallXea{L) 
^{z)-exp{z) = 0{zP+^). 

We solve Problem [T] using a convex optimization approach and bisection with respect to the 
CFL number z/, as described in [21] . This approach allows us to optimize methods with large 
numbers of stages in order to improve the maximum absolutely stable time step i^stab- The 
optimization is carried out for 2"'^- to 5**^-order accurate schemes; the constraint points A 
are taken as the spectrum of the SD semi-discretization of the same order of accuracy. 

Figure [s] shows the stability regions of the classical 4-stage 4*'^-order (ERK(4,4)) [26] and 
optimized 18-stage 4*'^-order (ERK(18,4)) methods superimposed on the scaled spectrum 
(i^stab K) of the 4*'^-order SD methods. Clearly, the optimized method allows the use of a 
much larger step size. Notice also that the stability region of the optimized method has 
nearly the same shape as the convex hull of the SD spectrum. 



3.2 Determination of Runge— Kutta coefficients 

The choice of stability polynomial does not fully determine the method; an ERK method of 
s stages has s(s + 1)/2 coefficients and only s of them are constrained by the stability polyno- 
mial. We now consider the problem of finding the RK coefficients A, b corresponding to a set 
of prescribed stability polynomial coefficients Pj. We use the remaining degrees of freedom 
to satisfy additional nonlinear order conditions, to obtain a low-storage implementation, and 
ensure that the truncation error coefficients are not too large. 

While the linear accuracy of the method is determined by the stability polynomial, the 
nonlinear accuracy depends on additional order conditions r- {A, b) = [5]. We ensure that 
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Re(z) Re(z) 

(a) ERK(4,4). (b) ERK(18,4). 

Figure 3: Examples of stability region for 4*'^-order ERK methods (red lines) and scaled spectrum 
(i^stab Aj for i = 1, . . . , A^DOJ") of the 4*''-order SD scheme (blue dots). 



those conditions are satisfied up to order p and we seek to minimize the Euchdean norm of 
the truncation error coefficients of order p + 1 [5l [16] : 

Memory requirements for RK methods are typically on the order of s x N^'^^ . To avoid 
the need for large amounts of memory, we employ the low-storage algorithm presented in [I23j , 
which reduces the requirement to 3 x A^^^^, i.e. 3 registers per stage. The coefficients of the 
methods are provided in terms of the low-storage formulation, which is given in Algorithm 
[TJ This algorithm also retains the previous solution values so that a step can be restarted if 
a prescribed stability or accuracy condition is not met. 

The optimization problem may be stated formally as follows 
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^3 ^ 

for 2 = 1 : s do 

t 4- t" + Q At 

5*2 ^ 5*2 + 5i 5*1 

Si ^ 71, 5i + 72, ^2 + 73, ^3 + A At 
end for 

^ Si 

Algorithm 1: Low storage implementation (3S*) 

Problem 2 (RK method optimization) 

maximize C*^^"*"^^ 
subject to 

Ti^-''^(A,b) = (0<j<p) 
b^A^-^e = /3j (0 < J < s) 
r(A,b) = 

Here r(A,b) = represents the conditions necessary for the method to be written in low- 
storage form. In practice, we impose those conditions implicitly by taking the low-storage 
coefficients as decision variables and computing the Butcher coefficients (A, b) from them. 

We use the RK-opt ("Runge-Kutta optimization") package to search for optimized meth- 
ods. This software uses MATLAB's f mincon function with the interior-point algorithm and 
the multi-start global optimization toolbox. Six hundred random initial guesses were used 
to find each optimized RK method. The RK-opt package and its extensions [2S] are freely 
available at ht tps : //github . com/ketch/RK-opt[ 

3.3 Efficiency and CFL number 

Time integration with an explicit method always incurs a step size restriction v < Vst&h 
based on stability. Accuracy typically also leads to a constraint on the time step, which for 
hyperbolic problems translates to a constraint on the CFL number, of the form u < z/acc, 
where face is the largest CFL number satisfying a prescribed error tolerance. Other concerns, 
such as positivity, may further restrict the CFL number, but we focus on z/gtab and z/acc- 

If '^stab < ^acc, Stability is the more restrictive concern and the relative efficiency of two 
RK methods of order p can be measured as the ratio of the maximum effective stable CFL 
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number p^tah/s: 

where a denotes a safety factor applied to both schemes. If Xstab > 1 then method 1 is more 
efficient. This quantity measures the relative time interval integrated per unit work [22]. 

On the other hand, if z/acc < t'stab, accuracy is the more restrictive concern so relative 
efficiency should be based on the ratio of step sizes, giving an equivalent global error. A first 
estimate, assuming that local errors simply accumulate, yields the relative efficiency measure 

1 

Xacc = I -TZ-TT (30) 



where C'f^^^ and C^^^^ are the principal error norms of the two RK schemes (see Section 



3.2 ). Note that this measure is meaningful only if both schemes have order p. Although (30 ) 
is probably too simplistic because the error at each time step feeds back into the computation 
at the next step, it is used as a guideline for the selection of RK schemes among the optimized 
methods presented in this papeij^ 

Figure |4] shows both Xstab and Xacc for 2^^^^- to 5*^-order optimized schemes over widely 
used traditional explicit ERK methods of the same accuracy: the mid-point rule ERK(2,2); 
Heun's 3-stage 3'''^-order ERK(3,3) method [H]; the classical 4-stage 4*'^-order ERK(4,4) [26j; 
and the 6-stage S^'^-order Runge-Kutta-Fehlberg ERKF(6,5) method [12]. The maximum 



stable CFL number i/gtab defined in Section 2.2.1 is also shown for completeness. 



3.4 Discussion 

We provide optimized methods for s < 20, because for larger values of s, the convex solvers 
used in the algorithm of [21] often fail due to poor numerical conditioning. However, Fig- 
ures 4(a) and 4(b) already show that the marginal efficiency gain achieved by adding another 
stage, becomes vanishingly small for large s. 

Indeed, the asymptotic efficiency gain that could be achieved by using additional stages 
is bounded, since the classical CFL theorem implies that the scheme cannot be stable for a 
CFL number greater than s [33]. An even tighter bound can be inferred by recalling that 
the stability region of an s-stage ERK cannot contain the closed disk with diameter [—2 s, 0] 
as a proper subset [21j. By determining the largest of such disks contained in the spectrum 



^Notice that (30) differs from Equation (28) in [52] by tlie exponent. In Kennedy et al. [12] tfie power of 
C2^^^^ / c[^^^^ is set to l/{p + 1) because the control of the local truncation error is the main objective of 
the comparison [18) . 
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(a) Optimal 2"'^-order methods vs. ERK(2,2). 
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(b) Optimal B'^'^-order methods vs. ERK(3,3). 
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(c) Optimal 4*''-ordcr methods vs. ERK(4,4). 
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(d) Optimal 5 -order methods vs. ERKF(6,5). 



Figure 4: Efficiencies and maximum linearly stable CFL number of tlie optimal ERK methods 
over some traditional ERK schemes of the same accuracy. 
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Method 


s 




Xstab Xacc 


Midpoint ERK(2,2) 


2 


1.7678x10^°^ 1.7180xlO~^^ 


1 1 


Optimal ERK(3,2) 
Optimal ERK(8,2) 


3 
8 


1.9587x10-°! 7.5938xl0-°2 
2.0968x10-°! 1.1294xl0-°2 


1.11 1.00 
1.19 0.98 



Table 1: Step size, stability efficiency Xstab find estimation of accuracy efficiency Xacc of the 
selected optimal 2'^'^-order ERK methods. The reference 2-stage method (corresponding to the 
values in italics) is the midpoint ERK method. 



Method 


s 




Xstab Xacc 


Heun's ERK(3,3) 


3 


7.5755xlO-°2 4.6296xl{]-^'^ 


1 1 


Optimal ERK(5,3) 
Optimal ERK(18,3) 


5 

17 


9.0719xl0-°2 9.9290x10-°^ 
1.0718x10"°! 7.1115x10"°^ 


1.20 1.00 
1.42 0.71 



Table 2: Step size, stability efficiency Xstab and estimation of global accuracy efficiency Xacc of 
the selected optimal S'^'^-order ERK methods. The reference method (corresponding to the values 
in italics) is Heun's 3-stage method. 

of the SD method, upper bounds on the efficiency of optimized methods can be obtained. 
A further refinement can be obtained by using Theorem 5 of [33], which refers to ellipses 
instead of only disks. 

These considerations imply that, for hyperbolic PDE discretizations, only a limited num- 
ber of stages are necessary to realize most of the potential efficiency gain. Based on Figure 
|4| it seems that the number of stages that provides a significant improvement increases with 
the order p. 

The blue lines indicate that the global error efficiency Xacc of our schemes generally 
decreases with the number of stages. However, for all schemes, Xacc is within 40% of the 
reference scheme value. We emphasize that accuracy is not the primary concern in the design 
of these schemes; certainly better accuracy could be obtained if one were willing to sacrifice 
some stability. 

Tables [T] to |4] list the value of Xstab and Xacc of the optimized schemes that are used for 
the test problems in the next section. Two methods have been selected for each order of 
accuracy. Those with fewer number of stages have an accuracy efficiency close to that of the 
reference methods, whereas the schemes with a large number of stages are characterized by 
a large value of the stability efficiency and an accuracy efficiency which is greater than 0.7. 
The coefficients of the selected optimized methods in terms of the low-storage formulation 
(see Algorithm [T] ) are listed in Appendix |Aj 

Since most of our new schemes use many stages, a thorough analysis of their internal 
stability [37| properties has also been performed. All the schemes are internally stable. 
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Method 


s 




Xstab Xacc 


Kutta's ERK(4,4) 


4 


3.9534x10^^'^ 1.4505x10-^^ 


1 1 


Optimal ERK(9,4) 
Optimal ERK(18,4) 


9 

18 


5.6977xl0-°2 5.0640x10"°^ 
6.5233xl0-°2 1.1087x10"°^ 


1.44 1.03 
1.65 0.75 



Table 3: Step size, stability efficiency Xstab and estimation of accuracy efficiency Xacc of the 
selected optimal 4*'^-order ERK methods. The reference 4-stage method (corresponding to the 
values in italics) is the classic ERK method. 



Method 


s 




Xstab Xacc 


Fehlberg ERK(6,5) 


6 


2.6916xl{)-^'^ 5.5557xlO-°3 


1 1 


Optimal ERK(10,5) 
Optimal ERK(20,5) 


10 
20 


3.6164xl0-°2 5.0975x10"°^ 
4.2195xl0-°2 1.0490x10"°^ 


1.34 1.39 
1.57 0.95 



Table 4: Step size, stability efficiency Xstab and estimation of accuracy efficiency Xacc of the 
selected optimal S^'^-order ERK methods. The reference method (corresponding to the values in 
italics) is Fehlberg's 6-stage, 5*^-order method. 

4 Applications 

In order to asses the efficiency and the accuracy of our new ERK schemes, we have performed 
a series of numerical simulations. The computations run on a machine with 2 x 2.4 GHz Quad- 
Core Intel Xeon, using the Coolfluid 3 collaborative simulation environment [35]. Sixteen 
gigabytes of RAM were available. The grids have been generated using Gmsh software jl4j . 
In Coolfluid 3, the CFL-number in two dimensions is defined as 

with ax and ay the x- and y-components of the wave speed a, and with Ax and Ay the width 
and height of a Cartesian grid cell. This definition, when applied to the Roe scheme on a 
structured Cartesian mesh leads to a positive and stable discretization for a CFL number 
smaller than unity [TU]. In practical computations it has been observed that slightly larger 
time steps may be used without affecting stability [TD]. 

4.1 Order verification 

In this section we present the convergence study of the optimized ERK scheme listed in 
Tables [T] to |4j We integrate a system of nonlinear non-autonomous system of ffist order 
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ODE [33122] 




(32) 



with the time t ranging from t° = 1 to t*^ = 1.4, and with the following initial conditions: 
= 1, (3'2(^°) = The analytical solution of this system is qi(t) = l/t, q2{t) = e~* . 
We use the norm of the error 



to study the time integration error. Here Qi and Q2 denote the numerical solutions. Figure 
[5] shows the norm of the error \e\ as a function of the time step At. It can be seen that for 
all ERK schemes the expected order of accuracy is achieved. Moreover, the new optimized 
methods show significantly smaller errors than the reference methods, as expected based on 
their smaller error constants. 

4.2 Advection of a Gaussian wave in an annulus 

The second problem we consider is the advection of a Gaussian wave in a 2D annulus. Such a 
problem models the transport of a scalar conserved variable q with variable advection speed 
u. The conserved variables and the convective flux are then 



Therefore, the conservation law reads 



where in a 2D Cartesian space 



u{x,y) = 



u 




-y 




= (jj 


V 




X 



Here u is the angular velocity which is set to = 27r. Note that the RK schemes have been 
optimized for the 2D advection equation with constant convective velocity whereas here a 
variable velocity is used. Therefore, the numerical results and performance presented in 



(gi(0-gi(0) + (g2(0-g2(0) 



dq 
di 



+ M . Vg = 0, 



(34) 
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this section can already be used to partially asses the robustness of the new time stepping 
methods. 

The initial Gaussian wave condition is centered at Xc = 7.5, i/c = 0.0 and is defined as 

qif) = e ^ , (35) 

where the radius is set to 6 = 0.6. 

The annulus is characterized by an internal radius = 5 and an external radius = 
10. One-quarter of the annulus is discretized for the actual computations. Simulations are 
performed using 2'^'^- to 5*^-order spatial and temporal discretizations from t° = to t*^ = 0.25 
(see Figure |6]). Several CFL numbers ranging from 0.05 to the maximum linearly stable one 
are used. A mesh with 110 x 110 (radial direction x azimuthal direction) quadrilateral cells 
with a maximum aspect ratio of 1.8 are used for the second-order computations. Such a mesh 
leads to a total number of DOFs is 48400 which is held constant for higher order accurate 
simulations by coarsening the grid in both directions. Extrapolation boundary conditions 
are imposed on both circular boundaries. 



10 

5 




-10 





I I 

10 -5 5 10 



0.5 







Figure 6: Gaussian wave advected in the annulus att = 0.25; solution computed with the 4*'^-order 
SD method and the optimal ERK(18,4) scheme. 



The exact solution is just a rotation of the initial solution and is given by (35) with 
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Xc = 0.0, = 7.5. 

We solve this problem with each of the reference and optimized ERK schemes using 
the predicted maximum stable CFL number for each scheme. In every case, the resulting 
computation is stable, confirming the theoretical prediction. Figure [7] shows the maximum 
norm of the error vector 

||£(r)||i^ =max|£,| = |max|(Q,(t'^)-ge:.,.(t'))| for z = 1, 2, . . . , iV°OF (36) 

and the CPU time for each scheme. Although the number of DOF is the same in all the 
simulations, the error decreases rapidly with increasing order of the discretization. Remark- 
ably, a unit increment of the order of accuracy leads to a reduction of the error of one order 
of magnitude and to a faster simulation (Figure [7(a)[ ). This shows the benefit of using high- 
order accurate methods for wave propagation problems. As predicted, some of the highly 
optimized ERK schemes yield somewhat larger errors. Figure |7(b) highlights the speed-up 



obtained with the optimized RK schemes over the standard methods for high-order accurate 
simulations. Indeed, for 4**^- and 5*'^-order computations the new schemes reduce the com- 
putational time by 40% and 38%, respectively. These values match very well the theoretical 
results shown in the previous sections. 

Figure [8] shows the maximum norm error as a function of the one-step effective CFL 
number. Here qex,i is the exact solution at the solution point (or DOF) i. 

Interestingly, for some methods the error increases with increasing CFL number, while 
for others it actually decreases. The latter behavior is reminiscent of the behavior of many 
low-order schemes that are more accurate for CFL numbers close to 1 and more dissipative 
for small CFL numbers. In addition, we point out that a combination of a quasi uniform 
grid and the maximum CFL numbers obtained during the first optimization step results in 
stable full discretizations. 



4.3 Acoustic wave propagation 

In this example we solve the linearized Euler equations (LEE), which model the propaga- 
tion of small perturbations in a mean flow field. They are frequently used to compute the 
propagation of acoustic waves in the absence of acoustic sources, e.g. turbulence production. 
They have been successfully used to solve in a hybrid approach for cavity flow [39], jet noise 
[1], and vortexblade interaction [6j. 

The LEE are derived from the compressible Euler equations which mathematically de- 
scribe the three physical conservation laws (i.e. conservation of mass, conservation of momen- 
tum and conservation of energy) for an inviscid fluid. Thus, the definitions of the conserved 
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Figure 7: Error and CPU time for the advection problem. The label s,p for eaeh point indicate 
the number of stages s and the order p of the corresponding scheme. Open circles are used for the 
reference methods; closed circles are used for the optimized methods. 



variables q and the flux vector f = [f g h] are 
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(37a) 



(37b) 



In these equations, p is the mass density, m, v and w are the Cartesian velocity components, 
p is thermodynamic pressure, and E is specific total energy. The velocity vector u\% \iiv w"^ 
and its magnitude is denoted by l^l. For an ideal gas, which approximates well the ther- 
modynamic behavior of air in a wide range of thermodynamic conditions, the specific total 
energy E is related to the pressure and the velocity field by 



E^ 



1 p u} ^v"^ + w^ 



+ 



(38) 
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Figure 8: Influence of the CFL number on the maximum error norm of a 2D Gaussian wave 
advecting in an annulus. 
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where 7 = 1.4 is the heat capacity ratio for air. Equation (38) closes the hyperbohc system 
(37) of five nonhnear PDEs with five unknowns. 

The LEE are obtained from (37) by decomposing the primitive flow variables p, u and p 
into a mean flow value (■)o and a perturbation about this mean flow (■)': 



P 
u 

P 



Po + P', 

Mo + u', 

Po+p'. 



(39) 



Substituting these relations in (37), subtracting the mean flow terms and neglecting 
products of perturbations, the following sets of conserved variables q and the flux components 



f = [f g h] are obtained 



v 



(40) 



/ P' 

Pou' 
Pov' 
Pow' 

p' 

( Pqv' + vqp' \ 
PoVqu' 

PqVqv' + p' 

PoVqw' 
V Vop' + 'jpov' J 

This procedure also leads to a source term involving mean flow gradients (right-hand side 

of (IB 



/ pqu' + uqp' \ 

PqUqu' + p' 

PoUqv' 

PqUqW' 
\ Uop' + 'jpou' J 



g 



/ pqw' + wop' \ 
PoWqu' 
PoWqv' 
PoWqw' + p' 
\ wop' + -fpow' j 



(41) 



/ \ 

(po«' + «op') '-^ + {p^V + t;op') ^ + (po^' + ^op') If 
(po«' + u^p') ^ + iPov' + vop') ^ + {poW + Wop') 
ipou' + uop') If + ipov' + vop') ^ + {pow' + wop') |f 



v 



(7 - 1) (l'' (t? + If + 1?) - - ^'t - ^'t) 



(42) 



which partially accounts for the refraction effects. The source term is zero in case of a 
uniform mean flow. 

The initial solution for this numerical test has a Gaussian profile centered at the origin 
of the axes and it is given by 

lO-^e 



P 
pi 
u' 



cl P' 



(43) 



0, 
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where the radius of the Gaussian pulse is set to b = 0.05. The uniform mean flow variables 
are po = 1, po = 1, •y = 1.4 and uq = vq = 0. Simulation are done from = to = 0.3. 
The numerical domain is a circle with radius r = 0.5, which is also centered at the origin of 
the axes. For the 2"'^-order calculations, a mesh with 180 x 45 (radial direction x azimuthal 
direction) quadrilateral cells with a maximum aspect ratio of 1.85 is used. Therefore the 
total number of DOFs is 32400. As for the previous numerical test, this number is kept 
constant for higher order accurate simulations by coarsening the grid in both directions. 
Simple extrapolation boundary conditions are used. Figure [9] shows the contour plot of the 
acoustic pressure field at t = = 0.3. 




-0.5 1 1 1 1 1 

-0.5 -0.25 0.25 0.5 



Figure 9: Acoustic pressure contour p' at t = 0.3; solution computed with the 4*^-order SD 
method and the optimal ERK(18,4) scheme. 

The exact solution for the acoustic pressure field p' = p — Po obtained by integrating the 
LEE is used as a reference solution to compute the numerical error. Its analytical expression 



is given by (44) 

p' (t, x,y) = ^ e-(^) cos (^cot) Jo (^r/) ^d^, (44) 
with 7] = \ [x — 0.5)^ + (y ^ 0.5)^ and Jq the Bessel function of the first kind of order zero. 
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Figure [10] shows the maximum norm of the error and the CPU time for each scheme, using 
the predicted maximum stable CFL number. All schemes are again stable at their respective 
theoretical CFL values. The results are similar to those shown in Figure [7] for advection, 
although it appears that the spatial errors are even more dominant for this problem as the 
overall error is nearly the same for all time-stepping schemes of a given order. 



Figure 11 shows the maximum norm of the error versus effective CFL number for a range 
of CFL numbers. We observe that for a fixed order of accuracy the error is almost indepen- 
dent of the CFL number. More precisely we find that as long as the time step is smaller than 
the theoretical maximum stable value, the error is dominated by the spatial discretization 
error. Note that a unit increment of the order of accuracy of the full discretization leads to 
a reduction of the error of one order of magnitude. The optimized RK schemes speed up the 
simulations considerably. 
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Figure 10: Error and CPU time for the acoustic pulse problem. The labels s,p for each point 
indicate the number of stages s and the order p of the corresponding scheme. Open circles are used 
for the reference methods; closed circles are used for the optimized methods. 



4.4 Vortex shedding past a wedge 

This test case focuses on the von Karman vortex street past a triangular wedge |10] com- 



puted with the compressible Euler equations (37). Indeed, in the inviscid framework, vortex 
shedding phenomena can be described when the considered bodies present sharp corners 
which ensure the separation of the flow. This numerical test represents a more realistic 
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Figure 11: Influence of the CFL number on the maximum error norm of a 2D acoustic pulse 
propagating in a circular domain. 
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Figure 12: Configuration of the wedge problem. 



application and it is used to study the performance of the new temporal schemes for a non- 
linear system of PDEs and highly unstructured mesh. The compressible Euler equations are 
generally used to model the flow of an inviscid fluid, or the flow of a viscous fluid in regions 
where the effects of viscosity and heat conduction are neghgible. Typical applications include 
preliminary aircraft design and rotor-flow computations. 



In Figure [12] the configuration of the test case is illustrated, where the incoming flow is 
from left to right. The wedge is placed on the centerline y = of the computational domain 
and it is characterized by a length L. At the left boundary (the inflow) the flow is prescribed 
to be uniform with zero angle of attack and free-stream Mach number of 0.2. Both inlet 
density and inlet pressure are set to one. A pressure outlet boundary condition is imposed 
on the right boundary of the domain which is placed about 15 L away from the wedge. 
Far-field boundary conditions (i.e. uniform Dirichlet boundary conditions for the conserved 
variables) are imposed both on the top and bottom boundaries. 

An unstructured grid with 11686 quadrilateral cells with a maximum aspect ratio of 
1.78 is used for the 2'^'^-order calculations. This leads to a total number of DOFs is 46744. 
The number of DOFs are again kept about constant for higher order accurate simulation by 
coarsening the grid. For this test the exact solution is not available. Therefore a reference 
solution is numerically computed by solving the problem on the mesh with 11686 quadrilat- 
eral cells with the 5*'^-order SD method (292150 DOFs) and the ERKF(6,5) scheme. A CFL 
number z/ = 0.1 is used for the reference computation. 

In order to avoid discontinuities near the surface of the wedge during the transitional 
phase that is produced by the uniform free-stream initial conditions, an intermediate solution, 
in which the formed vortices have not yet separated, is computed with I'^^-order SD and 
ERK(2,2). That solution is used as the initial condition for all higher-order computations 
(including the reference one) which are carried out from t° = to = 200 (see Figure 13) 
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to generate new unsteady laminar initial solutions for each order of accuracy. Afterwards, 





Figure 13: Density contour of the flow past a wedge at t = 200; solution computed with the 
4*'^-order SD method and the optimal ERK(18,4) scheme. 



starting with these intermediate solutions, several computations are performed using the CFL 
number z/gtab for each scheme and measuring the error after 0.1 seconds. Figure 14 shows the 
maximum norm of the error and the CPU time for each scheme. Remarkably, we observe 
that the new schemes, designed using linear advection on a uniform grid, perform very well 
for the compressible Euler equations on an unstructured grid. Indeed, they speed up the 
simulations considerably, while retaining a small error norm. Moreover, we highlight that the 
use of a quasi-uniform grid and the CFL number z/gtab results in stable full discretizations. 
Therefore, also for this nonlinear test the theoretical stability efficiencies obtained in the first 
optimization step by using the 2D linear advection equation model are recovered. 



5 Conclusions and future work 



In this work we have developed new robust optimized explicit Runge-Kutta schemes for the 
spectral difference method to efficiently and accurately solve wave propagation problems on 
unstructured uniform or quasi-uniform non-simplex cell grids. We have shown that by using 
low-storage schemes with optimized stability function and reasonable leading truncation er- 
ror constant, one can significantly improve the performance of the resulting method of lines 
discretization. By integrating high-order accurate spectral difference semi-discretizations 
(i.e. 3^'^-, 4**^- and S^'^-order) with optimized Runge-Kutta methods, we have found stability 
efficiency improvements of 42% to 65%, for typical systems of hyperbolic conservation laws 
used in fluid dynamics. These performance gains correspond to a reduction in computational 
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Figure 14: Error and CPU time for the wedge problem. The labels ,s,p for each point indicate 
the number of stages s and the order p of the corresponding scheme. Open circles are used for the 
reference methods; closed circles are used for the optimized methods. 



cost of 29% to 40% for a fixed simulation time. These improvements, which agree remark- 
ably well with theoretical predictions based on analysis of the 2D advection equation, are 
obtained without significant sacrifices in accuracy. Indeed, when both spatial and temporal 
discretizations have the same order of accuracy the spatial error typically dominates. There- 
fore, for a fixed order of accuracy the error is almost independent of the CFL number and 
large time steps can be used. 

Our results also highlight the advantage of high order schemes, which is even more pro- 
nounced when optimized time integrators are employed. The schemes designed in this work 
are intended for the solution of purely hyperbolic systems of conservation discretized on 
unstructured grids. The goal of ongoing research involves the optimization of such a family 
of schemes for convection-dominated problems with diffusion, where anisotropic grids are 
needed to discretize the domain with an economical distribution of cells. 

Finally, we want to highlight that by keeping the number of degrees of freedom constant, 
while increasing the order of accuracy, the new schemes allow to get much more accurate 
solutions in about the same time, with about the same memory requirements. For example, 
using our optimized schemes, the proposed S^'^-order discretization is actually faster than 
the 2°^-order discretization, when the number of degrees of freedom is held constant. Of 
course, the S^'^-order discretization is also much more accurate. We expect that similar 
improvements could be obtained by using our approach to design optimized ERK schemes 
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for other high-order accurate semi-discretizations. 

A Low storage ERK coefficients 
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2.3002859824852059X 10-°' 
4.0500453764839639 x 10-°' 
8.9478204142351003x10-"' 
7.2351146275625733x10-"' 


2.3002859824852059x10-"' 
3.0214498165167158x10-°' 

8.02,560102,388.56679x10-°' 
4.36216188715117,53x10-"' 
1.129270.5979513513x10-"' 


O.OOOOOOOOOOOOOOOOx 10+™ 
2.5876919610938998X 10-°' 

-1.3243708384977859x10-°' 
5.0556618948362981 X 10-"^ 
5.670.5,50788,3024708x10-"' 


1.0000000000000000x10+™ 
5.5284013909611196x10-°' 

6.731851,3326032769x10-°' 
2.80310-54965521607x10-"' 
6.521511.581,59187,58x10-"' 


0.0000000000000000x10+™ 
0.0000000000000000x10+™ 

0.0000000000000000x10+°° 
2 . 75257979463342 1 3 X 10-°' 
-8.9.50.544,502214,8511x10-°' 


1.0000000000000000x10+™ 
3.4076878915216791 x 10-°' 
3.4143871647890728x10-°' 
7,2292984084963262x10-°' 
O.OOOOOOOOOOOOOOOOx 10+™ 



Table 7: 3S* low storage coefficients of ERK(5,3) method, optimized for 3'''^-order SD scheme. 
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c 


,a 


''1 


72 




fS 


0.0000000000000000x10+"" 
4.9565403010221741 x lO""' 
1 .3068799001687578 x 10""! 

-1.5883063460310493X 10-"' 
3..56811447401!)0!)3.^)Xll)-"' 
7.6727123317642698x10-"^ 
1. 0812579255374613 X 10-°' 
1.8767228084815801 x 10"°' 
9.6162976936182631x10-"' 

-2.2760719807560897x11)-"' 
1.1115681606027146x10+"" 
6.1266846427676520x10-°' 
1 .0729473245077408 x 10+°° 
3.7824186468104548x10-"' 
7.!)041891347640720x 10-"' 

-1.0406955693161675X 10+"" 

-2.4607146824557105x10-°' 

Table 8: 3S* 1 

c 


4.9565403010221741x10-"' 
9.7408718698159397X 10-°' 
-1.7620737976801870x10-°' 
1.4852069175460250X 10-"' 
-3.3127657103714951x10 "' 
4.8294609330498492X 10-"' 
4.9622612199980112x10-°' 
8.7340766269850378x10-°' 
-2.8692804399085370x10-"' 
1.2679897532256112x10+"" 
-1.0217436118953449x10-"' 
8.4665570032698360x10-°' 
2. 8263854742688246 X 10-°' 
-9.2936733010804407X 10-"' 
-8.47!)81247008l):!512 X 10-"' 
-1.6923145636158564X 10-"' 
-4.7305106233879967X 10-°' 

3w storage coeffi 


0. 0000000000000000 X 10+"" 
7.9377023961829174X 10-°' 
-8.3475116244241754x10-°' 
-1.6706337980062214X 10-°' 
3.04100915003:11427x11) "' 
6.9178255181542780X 10-"' 
1 .48871 16004739182 x 10+°° 
4.5336125560871188x10-°' 
- 1 . 27057760464587.39 x 10-°' 
8.3749845457747696x11)-"' 
1.5709218393361746x10-"' 
-5.7768207086288348x10-°' 
-5.7340394122375393x10-°' 
- 1 . 2050734846614470 x 10+°° 
-2.8100719513641002x10+°° 
1 .6142798657609492 x 10-°' 
-2.5801264756641613x10+°° 

cients of ERK(1' 

7i 


1.0000000000000000x10+"" 
3.2867861940811260x10-°' 
1.1276843361180819x10+°° 

1.314944739.5238016x10+°° 
5.20628915:!4209055 X 10-"' 
8.8127462325164985X 10-"' 
4.2020606445856712x10-"' 
7.6532635739246124x10-"' 
4.4386734924685722X 10-"' 
6.0.5030!):l95519!)(iS2xll)-'" 
1.5850209163184039x10+"" 
1.1521721573462576x10+°° 
1 . 1 172750819374575 x 10+°° 
7.7630223917584007X 10-°' 
1.0046657000652295 x 10+"" 
-1.9795868964959054X 10-"' 
1 .3360583594705518 x 10+°° 

7,3) method, opt 

72 


0.0000000000000000x10+"" 
O.OOOOOOOOOOOOOOOOx 10+"" 
O.OOOOOOOOOOOOOOOOx 10+"" 
8.4034574578399479X 10-"' 
8.5047738439705145 x 10-"' 
1.4082448501410852X 10-"' 

-3.2678802469619369X 10-"' 
5.3716357620635535X 10-"' 
9.0228922115199051x10-"' 
1.5960226946983552 x 11)-"' 
1.1038153140686748x10+"" 
1 .0843516423068365 x 10-°' 
4.6212710442787724x10-°' 

-3.3448312125108,398x10-°' 
1. 1 153826507090690 x 10+°° 
1 . 5503248734613539 X 10+°° 

-1.2200246424704212x10+°° 

imized for S'^'^-oi 

73 


1 . 0000000000000000 X 10+"" 
-3.7235794357769936x10-°' 

3. 33 15440189685536 x 10-°' 
-8.2667630338402520x10-°' 
-5.4628377681035534x10-°' 

6. 02 10777634642887 X 10-°' 
-5.7528717894031067x10-°' 

5.0914861529202782x10-°' 

3.8258114767897194x10-°' 
-4.027906322118,5290x10-°' 
-2, 0820434288562648 X 10-°' 

1 .4398056081652713 x 10+°° 
-2.8056600927348762x10-°' 

2.2767189929551406x10+°° 
-5.8917.5:i()10(|.54(i356xlO-°' 

9.1328651048418164x10-°' 

0.0000000000000000x10+°° 

der SD scheme. 

s 


0,0000000000000000x10+°° 
2,8:i(i:!4:i24S1011769xlO-°' 
5.4840742446661772X 10-°' 
3.6872298094969475x10-°' 
-6.8061183026103166x10-°' 
3.518,526,585510.5619x 10-°' 
1.6659419385562171 x 10+°" 
9.7152778807463247x10-°' 
9.0515694340066964x10-°' 

Table 9: 3S* 


2.836343248101 1769 x 10-°' 
9.73649807474X(i4(i:!xl0-°' 
3.3823592364196498x10-"' 
-3.6849518935750763 x 10-°' 
-4.1139587669859462x10-'" 
1.4279689871485013x10+°° 
1.8084680619636603x10-°' 
1.6057708856060601 x 10-°' 
2.9522267863254809 x 10-°' 

ow storage coeff 


O.OOOOOOOOOOOOOOOOx 10+°° 
-4.655641:!8:i7.5Gl:!01 x 10+"" 
-7.7202649689034453X 10-"' 
-4.0244202720632174x10+°° 
-2.1296873883702272x10-°' 
-2.4350219407769953X 10+°° 

1.9856336960249132x10-'" 
-2.8107894116913812x10-°' 

1 .6894354373677900 X 10-°' 

icients of ERK(9 


1. 0000000000000000 X 10+°° 
2.4992627683300688 x 10+°° 
5.8668202764174726x10-°' 
1.2061419816240786x10+°° 
3.4747937498564541 x 10-°' 
1.3213458736302766x10+°° 
3. 1 1903(i;i45;i264964 x 10-°' 
4.3514189245414447X 10-°' 
2.3596980668341213x10-°' 

,4) method, opti 


0.0000000000000000 x 10+°° 
0.000000001)001)001)0x10+°° 
0.0000000000000000 x 10+°° 
7.6209857891449362x10-"' 
-1.9811817832965520x10-"' 
-6.2289587091629484 x 10-°' 
-3.7522475499003573x10-°' 
-3.3554373281046146 x 10-°' 
-4.5609629702116454X 10-"' 

mized for 4*^-or 


1 .0000000000000000 X 10+°° 
1.2029238731008208x10+°" 
7.5749675232391733x10-°' 
6. 1636907196196419X 10-°' 
-2.7463346616574083x10-°' 
-4.3826743572318672x10-°' 
1.2735870231839268X 10+°" 
-6.2947382217730230x10-°' 
0.0000000000000000x10+°° 

der SD scheme. 


O.OOOOOOOOOOOOOOOOxK)!"" 
1.238416!)48(Ki202!)8xl() "' 
1.1574324659554065X 10+"" 
5.4372099141646926x10-°' 
8.8394666834280744x10-°' 
-1.221204217660,5774x10-"' 
4.4125685133082082 x IQ-"' 
3.8039092095473748x10-"' 
5.4591107347628367x10-°' 
4.8731855535356028x10-°' 
-2.3007964303896034x10-"' 
-1.8!)()7()5(i(i()291.5873xl0-"' 
8.1059805668623763x10-"' 
7.7080876997868803x10-°' 
1.1712158507200179x10+°° 
1.275535101800.1.545 x 10+"" 
8.0422.507946108.564x10 "' 
9.7.5(180802.50761848x11)-"' 

Table 10: 3S* 

c 


1.2:i841(i!)48062029«xl() 
1.01762025.14280349x11)+"" 
-6.9732026387527429X 10-"' 
3.4239356067806476x10-°' 
1 .8177707207807942 x 10-°' 
-6.1188746289480445X 10-°-'' 
7.8242308902580354x10-"' 
-3. 764286475053295 1x10-"' 
-4.5078383666690268x10-°' 
-7.5734228201432585x10-°' 
-2.7149222760935121x10-°' 
1.1833084341657344X1I)-"' 
2.8858319979308041 x 10-"' 
4.6005267586974667X 10-°' 
1 . 8014887068775631 X 10-°' 
-1., 5,50817539,5461857x10-°' 
-4,00!)5737929274988x 11)-"' 
1,4949078307038011x10-"' 

OW storage coeff 


O.OOUOOOUOOOOOOOOOxK)!"" 
1.175081!)Hll!)51678xll)t"" 
3.0909017892654811x10-"' 
1.4409117788115862x10+°° 

-4.3563049445694069x10-°' 
2.0341503014683893 x 10-°' 
4.9828356971917692x10-°' 
3.5307737157745489x10+°° 

-7.9318790975894626 x 10-°' 
8.9120513355345166x10-°' 
5.7091009196320974x10-°' 
1.691218857,501.5419x10-°' 
1.0077912519329719x10+°° 

-6.8532953752099512x10-°' 
1.0488165551884063x10+°° 
8. 364776 1 371829943 X 10-°' 
1.3087909830445710x10+°° 
9.0419081700177323x10-°' 

icients of ERK(1 

71 


1.00l)()001)(J(l01)()(l01)()xl()i"" 

-1.2891068509748144x11)-"' 
3.5609406666728954X 10-"' 

-4.0648076226104241 x 10-°' 
6.0714786995207426x10-°' 
1.0253501 186236846 x 10+°° 
2.4411240760769423x10-°' 

-1.2813606970134104x10+°° 
8.1625711892373898x10-°' 
1.0171269354643386x10-°' 
1.9379378602711269x10-°' 
7.440H043544851782X 10-°' 

-1.2591764563430008x10-°' 
1.1996463179654226x10+°° 
4.5772068865370406x10-°' 
8.3622292077033844X 10-°' 

-1.4179124272450148x11)+"" 
1.366145900.5331649x10-"' 

8,4) method, op 

72 


O.OOOOOOOOOOOOOOOOx 11)1"" 
0.00(100000000000(10x10+°° 
0. 0000000000000000 X 10+°° 
2.5583378537249163x10-°' 
5.2676794366988289x10-°' 

-2.5648376621792202x10-°' 
3.1932438003236391 x 10-°' 

-3.1106816010852862x10-°' 
4.7631196164025996x10-°' 

-9.8863727938895783x10-°' 
1.9274726276883622x10-°' 
3.23898608.5.5971.508x10-°' 
7.5923980038397509x10-°' 
2.0635456088664017X 10-°' 

-8.9741032556032867x10-°' 
2.68999,T25O5676190xl0-°' 
4.1882069379552307x11)-'" 
6.2016148912381761x10-'" 

timized for 4*^^-0 

73 


l.OOOOOOOOOOOOOOOOxlO"'" 
3.,581(i,500441970289xl0-"' 
5,8208024465093577x10-°' 
-2.2615286894283538x10-°' 
-2.1715466578266213x10-°' 
-4.6990441450888265X 10-°' 
-2.7986911594744995x10-°' 
9.8513926355272197x10-°' 
-1.1899324232814899x10-°' 
4.2821073124370562x10-°' 
-8.2196355299900403x10-°' 
5.811399705767,5074x10-°' 
-6.1283024325436919x10-°' 
5.6800136190634064X 10-°' 
-3.3874970570335106x10-°' 
-7.30712,38125137772x10-°' 
8,3936016900374532x10-°' 
0,0000000000000000x 10+°° 

rder SD scheme. 

s 


0.0000000000000000x10+™ 
2.5978836767039448x10-°' 
9.9045731168085657x10-'" 
2.16551 18823045644 x 10-"' 
5.0079500784155040x10-'" 
5.5922519148647800x10-'" 
5.4499869734044426x10-°' 
7.6152246625852738x10-°' 
8.4270620830633836x10-°' 
9.1622098071770008x10-°' 


2.5978835757039448x10-"' 
1.7770088002098183x10-°' 
2.48l(i:i0637:!161344xl0-°' 
7.9417368275785671 x 10-"' 
3.8853912968701337x10-"' 
1.4650516642704694x10-"' 
l.,587517,379465,5811xl0-°' 
1.6506056315937651 x U)-'" 
2.1180932999328042x10-°' 
1.5693923403495016x10-"' 


0.0000000000000000x10+"" 

4,043660078,5287713x 10-°' 
-8„50:!427404129,5027xl0-'" 
-6.9508941671218478x10+°° 

9.2387662262320684x10-"' 
-2.5631780399589106x10+"" 

2,.")457,i48n99988827xl0-°' 

3,125^,!17,l:!07014,54xlO-°' 
-7.0071148003175443x10-°' 

4.8396209710057070x10-"' 


1 .0000000000000000 X 10+"" 

6.8714670697294733X 10-°' 
1.0930247004585732x10-°° 
3.2259753823377983x10+°° 
1.0411637008416110x10+"" 
1.2928214888638039x10+"" 
7.3914627692888835X 10-°' 
1.2391292570651462x10-°' 
1.8427534793568445 x 10-°' 
5.7127889427161162x10-"' 


0.0000000000000000x10+"" 

O.OOOOOOOOOOOOOOOOx 10+°° 
O.OOOOOOOOOOOOOOOOx 10-°° 
-2.3934051593398129x10+°° 
-1.9028644220991284x10+"" 
-2.8200422105835639x10+"" 
-1.8,326984641282289x10^°° 
-2.1990945108072310x10-°' 
-4.0824306603783045x10-°' 
-1 .377669791 1236280 x 10-"' 


1 .0000000000000000 X 10+"° 

-l.,33177840914003,36xl0-°' 
8,2604227852898304x10-"' 
1.5137004305165804X 10+°° 

-1.3068100631721906x10+°° 
3.0366787893355149x10+°" 

-1.44945S2670S31953xl0+"° 
3.83431;!87.;3(,S5103xlO+"° 
4.1222939718018692X 10+°° 
0.0000000000000000x10+°° 



Table 11: 3S* low storage coefficients of ERK(10,5) method, optimized for 5* -order SD scheme. 
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O.OOOOOOOOOOOOOOOOx 
1.7342385375780556X 

3.0484982420032158X 

5.52713!)5fi4.")729193x 

4.7079204549750037 

1.5652540451324129X 

1.8602224049074517X 

2,8426620035751449x 

9,5094727548792208x 

6.8046501070096010X 

5.9705366562360063X 

1.8970821645077286X 

2.9742004004529606 

6.0Sl:l4(i:l700l:14940 

7.3080004188477765X 

9.1656999044951792X 

1.4309687554614530X 

4.1043824968249148 

8.4898255952298962 

3.3543890258348421 



To™" 
10-°' 

10-°' 
10-°' 
10-°^ 
10-°' 
10-°' 
10-°' 

II) 

10-°' 
10-°' 
10+°° 

10-°' 
10-°' 
10-°' 
10-°' 
10+°° 
10-°' 
10-°' 
10-°' 



10=01- 
10-01 

10-°' 
10-°' 
10-01 

10-02 

10-01 

10+00 

io-°2 

10-01 
10+00 

10-03 

10-°' 

10-02 

10-°' 
10+00 

10-04 

10-02 

10-°= 

10-02 



71 



lO+nr 
10+00 

10+00 

10-°' 

10-03 

10+00 
10+00 

10+°° 
10-°' 
10+°° 
10-01 
10-01 
10-°' 
10-°' 
10+°° 
10-01 
10-03 
10-01 

10-°' 
10+00 



72 



0000000000000000 X 10+"" 
8952052154583572x10-°' 

8988129100385194x10-°' 
5701564494077057x10-°' 
4232462479216824X 10-°' 
2727083024258155x10+°° 
1126977210342681 x 10+°° 
1360709645409097X 10-°' 
1181089682044856x10-°' 
7881272382085232x10-°' 
9032886260666715 x 10-°^ 
1871051065897870x10-°^ 
4002463790fl8e219x 10-°^ 
48972712511,547,50x 10-°' 
6244269699436817x10-°' 
7486056986590294X 10-°* 
3219312682036197X 10-°^ 
2852588972458059X 10-°' 
4473719351268962x10-°' 
4345446089014514X 10-°' 



0.0000000000000000 X 10+"" 
O.OOOOOOOOOOOOOOOOx 10+°° 

0.0000000000000000 X 10+°° 
1 .9595487007932735 x 10-°' 
^6.9871675039100595 x 10-°= 
1.0592231169810050x10-°' 
1.0730426871909635x10+°° 
8.9257826744389124x10-°' 
-1.4078912484894415X 10-°' 
^2.6869890558434262 x 10-°' 
-6.5175753568318007x10-°^ 
4.9177812903108553X 10-°' 
4.6017684776493678x10-°' 
-6.4689312947008251 X 10-°'' 
4.4034728024115377X 10-°' 
6.1086885767527943x10-°' 
5.0546454457410162 x 10-°' 
5.4668.509293072887X 10-°' 
7.14141S242l)99,5431xl0-°' 
-1.0,55809.5282893749x10+°° 



1x10+"" 
ixlO+»» 

10+°° 
10-°' 
10-01 

xlO-"' 
xlO-"' 

10-°' 
lO+oo 

10+00 

;xlo+»» 
X 10+00 

10-°' 
10+°° 
10+00 

.XlO-"' 
lxlO+»» 
10-01 

10-°' 
10+00 



1.7342385375780556X 
2.8569004728564801 X 

6.8727044379779589 

1.2812121060977319 

4.9137180740403122 

4.7033584446956857X 

4.4539998128170821 X 

1.2239824887343720 

2.0B16403986024421 

1.5941162575324802 

1.2963803678226099X 

1.7287352967302603X 

1.1660483420536467 

7.7997036621815521 

3.2563250234418012 X 

1.0611520488333197X 

6.5891625628040993 X 

8.3534647700054046 

9.8972579458252483 

4.3010116145097040 



O.OOOOOOOOOOOOOOOOx 
-1.1682479703229380X 

-2.5112155037089772X 
-5.3259900134735988X 

2.9243033509511740X 
-4.7948973385386493X 
-5.3095533497183016X 
-2.3624194456630736 

2.0008993730589547 
-1.4985808661597710X 

4.8941228502377687X 
-1.0387512755259576X 
-l.,328766427,3288191 

7.3858078822837511 
-4.3321586294096939X 

4.8199700138402146X 
-7.0924756614960671 X 
-8.8422252029506054 X 
-8.9129367099.345231 X 

l.,32971571,34040762x 



1.0000000000000000 
1.4375468781258596 

1.5081653637261594 
-1.4575347066002088 
3.1495761082838158 
3.5505919368536931 
2.3616389374566960 
1.0267488547302035 
3.5991243524519438 
1.5172890003890782 
1.8171662741779953 
2.8762263521436831 
4.6350154228218754 
1.3573122110727220 
2.0001066778080254 
9.1690694855534305 
2.0474618401365854 
-3.2336329115436924 
3.2899000734742177 
0.0000000000000000 



Table 12: 3S* low storage coefficients of ERK(20,5) method, optimized for 5*^-order SD scheme. 
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